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(N Abstract 

Although fully elastic when static, granular media become transiently elastic when being slowly 
sheared - during which both the elastic energy and stress relax. Starting from this observation, 
we cogently derive the framework for granular hydrodynamics, a set of differential equations con- 
c+h sistent with general principles of physics, especially reversible and irreversible thermodynamics. In 

o 

^ addition, an expression for the granular elastic energy is reviewed and further discussed. 
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I. INTRODUCTION 



In granular media, although the grains roll and slide, in addition to being compressed and 
sheared, only the latter, the deformation of the grains, leads to reversible energy storage 
that sustains a static, elastic stress, while rolling and sliding heats up the system. The 
granular strain field e^, therefore, has two contributions, an elastic one Uij accounting for 
deformation of the grains, and a plastic one for the rest, where = + p^. The 
elastic energy w{uij) is a function of u^, not of e^, and the elastic contribution to the stress 
Oij is given as ^-(ity) = —dw/duij. With Oij = tz^ in statics, stress balance Vy<7y = 
may be closed with = 7Tjj(wy) and uniquely determined employing appropriate boundary 
conditions [U [2] . Because the plastic part of the strain needed for arriving at a given stress 
state is quite irrelevant for its determination, one may with certain justification consider 
static granular media, say a sand pile at rest, as fully elastic. 

If this sand pile is perturbed by periodic tapping, circumstances change qualitatively. Its 
conic form will then degrade until the surface becomes flat. This is because part of the 
grains in the pile loose contact with each other temporarily, during which their deformation 
decreases. This implies a relaxing elastic strain u^, and correspondingly, smaller elastic 
energy w(iiij) and static stress Tiijiuij). Since the sand pile is no longer able to sustain static 
stresses, it is now a transiently elastic system, same as polymers - though the respective 
microscopic mechanisms are of course very different: temporary unjamming and rearrange- 
ment of the grains versus slow disentanglement of polymer strands. Note that flattening a 
sand pile implies sizable granular rearrangement, requiring a considerable portion of plastic 
strain p^. 

Quantifying the random motion of the grains as granular temperature T g , we may take 
the relaxation time r of the elastic strain function of T g , with r(T g ) — > oo for T g — > 0. 

For vanishing T g , there is no strain relaxation, the deformation of the grains are maintained, 
the sand pile keeps its conic shape, and the system is elastic. For finite T g , the elasticity 
turns transient, with u^, irij(v,ij) and w(v,ij) relaxing. 

When granular media are being slowly sheared, circumstances are similar. In addition to 
moving with the large scale velocity Vi, the grains also move and slip in deviation of it. This 
allows temporary, partial unjamming, and implies a finite T g , both again lead to transient 
elasticity. Since T g is not always an externally imposed parameter, as with tapping, but 
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frequently internally produced, especially by shear flows, it is an independent variable of the 
granular hydrodynamic theory, to be accounted for by its own equation of motion. More 
specifically, the production of T g by shear flows should have great similarities to viscous heat 
production in normal fluids. 

Granular media has different phases that, in dependence of the grain's ratio of deformation 
to kinetic energy, may loosely be referred to as gaseous, liquid and solid. Moving fast and 
being free most of the time, the grains in the gaseous phase have much kinetic, but next 
to none elastic, energy [3J. In the denser liquid phase, say in chute flows, there is less 
kinetic energy, more deformation, and a rich rheology that has been scrutinized recently [I] . 
In granular statics, with the grains deformed but stationary, the energy is all elastic. This 
state is legitimately referred to as solid because static shear stresses are sustained. If granular 
solid is slowly sheared, the predominant part of the energy remains elastic. As discussed, 
the system is transiently elastic, or quasi-solid. In this paper, we focus on the last two cases, 
and for simplicity refer to both as the solid granular phase. 

The transition between permanent and transient elasticity is a crucial key to understand- 
ing granular solids. And remarkably, it is as input quite sufficient for a formal and cogent 
derivation of the framework for granular solid hydrodynamics - if one takes careful notice 
of all general principles of physics, especially symmetry and thermodynamic considerations. 
This is the first part of the present paper. The second part deals with an concrete expression 
for the granular elastic energy, how this expression is supported by extensive experimental 
data from granular statics. This is important because general principles only confines the 
structure of the hydrodynamic theory - they yield a framework into which many different 
theories fit. The three sets of differential equations given below need the input of spe- 
cific expressions for the thermodynamic energy and the transport coefficients. Only when 
their functional dependence on the thermodynamic variables is given, do the theories attain 
predictive power. 

In the following, we first recall the hydrodynamic theory of permanent and transient 



elasticity, in § II A and § II B then merge both to form granular hydrodynamics, in § II C 



All equations in these three subsections are valid irrespective what form the energy w has. 
A specific energy density suitable for granular media is then reviewed and further discussed 
in 



III 



In an accompanying paper [5], we compare hypoplasticity [B], a state of the art engineering 
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model on the behavior of granular solids, with granular solid hydrodynamics as derived here, 
and specified using the elastic energy of § [TTT 



II. ELASTICITY THEORY 
A. Permanent Elasticity 

The conserved, thermodynamic energy density w of solids is a function of the symmetric 
strain field tty = Uji, and of the densities of entropy s, mass p, momentum g. So we write 
(neglecting gravity) 

dw = Tds + pdp + Vidgi — iiijduij, (1) 

denoting T(s,p,Uij) = dw/ds, p(s, p, g^Uij) = dw/dp, Vi = dw/dgi = gi/p, 7Ty(s, p,Uij) = 
—dw/duij. The equations of motion for the energy and its variables are 

§- t w + V.Q, = 0, f s + Vi/i = R/T, (2) 

y + Vdi = 0, §- t 9i + Vpij = 0, (3) 
f t Uij - v i: j + [\ViVj + u ik VjV k +i++j] = 0, (4) 

where 4 = + ffcVfc and = |(ViUj + VjUj). Expressing conservation laws and entropy 
production, the first four equations are quite general and shared by all hydrodynamic the- 
ories. Alone, they describe normal fluids and represent the simplest hydrodynamic theory. 
The fifth equation is characteristic of elastic systems, especially ones that break the trans- 
lational symmetry spontaneously. (More on why Eq Q must have the above form is given 
in [7].) Inserting Eqs (|2]|4]) into the temporal derivative of Eq Q, 

§t w = T M s + VmP + Vi-^gi-TTij^Uij, (5) 

and introducing the notations: ff, a®, taking them to be given as 

fi EE SVi ~ /f, (6) 

+(Ts + v i g i + fip + g i v j -w)-<Ty, (7) 

we obtain 

ViQi = Vi(Tfi + pji + VjOij - yjTTij) (8) 
+/f V;T + afjvij + ViVjitij - R. 
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Clearly, one can write the left hand side of Eq ^ as the divergence of something, plus 
something else that vanishes in equilibrium (because the so-called thermodynamic forces, 
VjT, Vij and Vj7Tjj do). Therefore, an inviting possibility is to identify the first with the 
energy flux Q iy and the second with the entropy production R, a quantity that also vanishes 
in equilibrium, 



Qi = Tft + nfc + VjOij - yj-Kij, (9) 
R = f t D ViT + a? Vij + //;V ; -,, (10) 

This identification is in fact unique. It is easy to verify that, as long as the energy w remains 
general, unspecified, there is no other way to write the left hand side of Eq (pjl) as the sum 
of a divergence and an expression that vanishes in equilibrium. 

Taking in Eq ( 10 ) (V«T, V{j, Vj7Ty) as the thermodynamic forces, (ff, a?, jji) as the fluxes, 
and forming each into a 12-component vector, Y and Z, the Onsager force-flux relation gives 
their linear connection as, 

Z = c-Y, (11) 

where c is the transport matrix, with diagonal elements that are positive, and off-diagonal 
ones that satisfy the Onsager reciprocity relation. The simplest example for c has only 
diagonal elements, all positive scalars, 

SF = «ViT, (12) 
<?$ = (vee$ij + VVij, (13) 
Vi = /3 P V i 7r ii . (14) 

Accounting for heat conduction and viscous stress, the first two equations are shared by all 
hydrodynamic theory. (The superscript °, here and below, denotes the traceless part of a 
tensor, eg. = % — ^SijVu.) The third accounts for permeation and defect motion, and is 



specific to elastic media [8], see section II A 1 



All elements of the matrix c, usually referred to as transport coefficients, are functions of 
the thermodynamic variables, s,p,Uij, or alternatively, of the conjugate variables, T, //, 7Ty. 
In the generally accepted and above employed linear version of the Onsager relation, they 
do not depend on thermodynamic forces, VjT, i>y, V_j7Ty. So we may take the coefficients 
k, rj, ( and f3 p to depend on the temperature, the pressure, and scalar combinations of the 
stress, such as 7r« and 7i 2 s = ^tt^-- 



1. Solid Creep Motion 



Enforcing a steady velocity at the surface of granular solid, the velocity field is observed to 
penetrate rather deep into the bulk of the granular medium, with a magnitude that decays 
exponentially with depth [9]. The usual collective modes of velocities in hydrodynamic 
theories of elastic media are of course such that they reduce to a constant velocity when 
stationary (sound), or one that varies linearly in space (shear diffusion). But there is also 



a less-known one that decays exponentially, a consequence of Eq ( 14 ) and the less studied 
permeation coefficient f3 p . We shall refer to this mode as "solid creep motion." 

Linearized with respect to velocity, Eq Q reduces, for the stationary case duij/dt = 0, 

to 

Vij = |/5 P Vfc(Vi7T jfc + VjlTik), (15) 

implying that mass and shear flows are possible without any changes in the elastic strain field, 
or equivalently, in the elastic stress and elastic energy. Similarly, momentum conservation, 
or Eq (|3]), linearized and in steady flow, d(pvi)/dt = 0, reduces to 

VjiDSy + Vij-W^) = (16) 

(where DSij stands for the diagonal terms that do not concern us here). Now, consider a half 
space y > filled with solid, which has its surface at y = 0, moving with a given velocity 
along x. Permitting only a ^/-dependence in this one-dimensional geometry, we have 



v xy = IfV^y, V y (7r xy - r,vl y ) = 0. (17) 



These two equations clearly imply exponentially decaying velocity v x and change of the 
elastic stress Sir xy , 

v x , Sn xy ~ exp - = . (18) 

In granular medium, this behavior will be modified, because the elasticity there is tran- 
sient rather than permanent. But should solid creep motion retains its qualitative behavior 



under certain circumstances, Eq (18) would constitute a natural explanation of granular 
creep flow. 
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B. Transient Elasticity 



Although the equations of the last section are fairly general and account for all kinds of 
elasticity, linear as well as nonlinear, they do exclude transient elasticity, such as realized 
in polymers. In these, + elasticity arises from entanglement of polymer strands, which are 
stretched and sheared, if not given enough time to disentangle. But if given enough time, 
the deformation, with it also the associated energy and stress, relax. So the system is to be 
accounted for by a set of equations which reduce to those of the last section for small time 
spans, but allow the deformation Uij to relax for longer time spans. 

The independent variables remain the same, so do the conservation laws. So Eqs ( Tpp ) 
are unchanged, but Eq Q is modified to allow for a relaxation term Xij 



at 



- Vij + [\ViVj + u ik VjV k + i j] = X i:j . (19) 



The same calculation of Eq (J5]), with the same notation of Eq (6]7), then leads to the same 
energy flux Qi, but a modified entropy production, 

R = if V 4 T + af> Vij + y&fKij + . (20) 

This implies 7Tj,- is now not only a conjugate variable, but also a thermodynamic force, 



increasing the dimension of the 12-component vector Y in Eq (11 ) by another 6 components. 



Similarly, the vector Z is also increased by the 6 components of X^, and c is now a 18 x 18 



matrix. Other from that, Eq (11) still holds. The simplest, diagonal and scalar example is 



again given by Eqs (12, 13, 14), in addition to 



Xij = pnl + p^eeSij, (21) 

a term that permits ity to relax, as long as is nonzero. 

As discussed in the last paragraph of § |II A the transport coefficients (3,(3\ are functions 



of the thermodynamic variable u^, or equivalently, of 7Tjj = TTij(uij). This remains true even 
though iTij is now also part of R, Eq (20), and hence an additional thermodynamic force. 



A point worth clarifying concerns the plastic strain p^: The total strain = +Pij, 
a purely kinematic quantity, obeys the equation 

f t Eij + [e ik VjV k + j]= v^. (22) 
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So, as a result of Eq (19), the plastic strain is determined by 



J] = 



(23) 



If a transiently elastic medium is quickly and uniformly deformed, such that there is no 
time for relaxation, Xij « 0, we have = Uij after the deformation. Holding it for a 
while, Vij = 0, the elastic deformation Uij relaxes to zero, while the plastic one pij grows 
accordingly, until one replaces the other completely, and we have = Eij. The system now 
stays where it is, and the initial displacement is referred to as "plastic," rather than elastic, 
because it does not have the tendency to return to the original position. 

Essentially this set of equations, as specified in this section, was recently shown well 
able of accounting for the full range of polymers' non-Newtonian behavior, including shear- 
thinning, elongational strain-hardening, rod climbing (the Weissenberg effect), and various 
empirical rules such as Cox-Merz and First Gleissle Mirror Rule [7J [10] 

C. Granular Elasticity 

As discussed, sand and other granular media display both elastic and transiently elastic 
behavior - depending on whether the granular temperature T g vanishes or not. Including 
the density of granular entropy s g as an additional, independent thermodynamic variables, 
the Gibbs relation of Eq (FT]) now reads 



Granular temperature is not a new concept. Haff, also Jenkin and Savage [3J, were 
probably the first to introduce it in the context of granular gas, using it to denote the 
average kinetic energy of the grains. Hence T g ~ e&, where is the kinetic energy density. 
Nowadays, this T g is routinely used in considering granular gas and liquid [IT]. Note that 
given this interpretation of T g , we have T g = dek/ds g ~ dT g /ds g , and the granular entropy is 
uniquely determined, s g ~ hiT s . More recently, there is much discussion of a configurational 
entropy S c in the literature. The original concept by Edwards was to approximate grains 
as infinitely rigid and all configurations as having identical energy [12], so S c is a function 
only of the system's volume. When relaxing the rigidity approximation, and allowing the 
elastic energy to vary, S c is again a function of energy and volume S c = S C (E,V), and a 
configurational temperature is naturally given as T c _1 = dS c /dE (see [13] for a review). 



dw = Tds + T g ds g + jjdp + Vidgi — TTijduij. 



(24) 
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In thermodynamics, the energy change dw from all microscopic, implicit variables is 
subsumed as Tds, with s the entropy and T = dw/ds its conjugate variable. From this, 
we divide out the kinetic energy of granular random motion, executed by the grains in 
deviation from the ordered, large-scale motion, denoting it as T g ds g , and calling s g and 
T g = dw/dsg granular entropy and temperature, respectively. In other words, we consider 
two heat reservoirs, the first containing the energy of granular random motion, the second 
the rest of all microscopic degrees of freedom, especially phonons. In equilibrium, T g = T, 
and s g is part of s. (In fact, we may simply forget s g , since it has far less degrees of 
freedom.) But when the granular system is being tapped or sheared, and T g is many orders 
of magnitude larger than T, then this leaky, intermediary heat reservoir can no longer be 
ignored. As s g then serves as a nonhydrodynamic, macroscopically slow degree of freedom, 
with T g its conjugate variable. 

Taking s g as the part of the entropy accounting for the granular kinetic energy, our 
definition is fairly close to the entropy of granular gas discussed above, as given by Haff, 
though its functional dependence will probably be modified, because it must be evaluated 
taking into consideration the effect of excluded volumes - an overwhelming one in the dense 
solid phase [see Eq (9) of the third of [llj]. The concept of configurational entropy, on the 
other hand, is closer to our second heat reservoir, the true entropy s, see section 6 of the 
first, and section 10 of the third, reference [2], for a discussion of their relationship. 

The functional dependence of s g (T g ), more precisely, the equation of state T g = 
T g (s, s g , p,Uij), is given once the energy w is known. Although all equations of this and 
the last two sections remain valid irrespective of what special expression one chooses for w, 
concrete predictions certainly depend on it. Since it appears difficult, at least at present, to 
evaluate w microscopically, one may alternatively employ experimental data in conjunction 
with general considerations to narrow down its possibility. We shall examine u>'s dependence 
on Uij and p in the next section, but defer that on s g to a future publication. 

Taking the balance equation for s g , in the uniform Sg = Rg/T g , we first of all 

need R g to contain the term — 7(T g — T) 2 . This is because being a slow, nonhydrodynamic 
variable, the equation of motion for s g should have the usual relaxation form, J^s 9 = — 7(T 5 — 
T), pushing T g towards the ambient temperature T. (Since any random motion of the grains 
implies such improbably high T g , neglecting T in this expression is always an excellent 
approximation. We shall therefore from here on always write §~ t s g = —^T g .) 
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Second, with the heat bath divided into two parts, viscous heat production should fill both 
baths simultaneously. Therefore, we keep the term cr^Vij in R, with a? = r/v^ + (vu, and 
write the analogous one, £^t> y into R g , with = rj g v^ + ( g vu denoting the viscous stress 
contribution from exciting granular random motion. The magnitude of the four viscosities 
depend on microscopic details and cannot be decided on general principles. For instance, 
while r] is probably a small quantity compared to rj g in dry sand, because macroscopic shear 
flows excite granular random motion first, rj should be quite a bit larger in wet sand: A 
macroscopic shear flow implies much stronger microscopic shear flows in the fluid layers 
between grains, and the energy dissipated in these layers should predominantly go to s, 
rather than to s g first. 

Third, granular entropy production R g should have the term K g VjT 2 , from an inho- 
mogeneous granular temperature, in exact analogy to the term fcViT 2 in R. So the final 
expression should be R g = Ti^Vij + K g VjT 2 — 7T 2 . A direct and desirable consequence of 
this expression is that for stationarity, ^s g = R g /T g = 0, and a constant T g , any shear 
flows excite the granular temperature of 7T 2 = rjgV^v^ + Cg v u, which is (as discussed) what 
renders granular elasticity transient. 

We do not have good reasons for ruling out a term in R g analogous to yiVjiTij, or one 
~ ViTg in R. But neither is there any experimental evidence demanding their existence. 
So although both are allowed for the general case, they are left out here for the simplicity 
of display. On the other hand, a term in R g analogous to Xijir^ cannot exist, because we 
would then have 7T 2 = X^ir^ for granular statics, implying a finite T g and decaying sand 
piles. 

Given the above consideration specifying R g , we may embark on the derivation of the 
equations of motion for granular elasticity, in the same way as above. We start from the 
following equations, 

f t w + ViQi = 0, y + Vdi = 0, (25) 
f s + V,/, = R/T, §- t s g + V t F t = R g /T g , (26) 

d_ 



dt 



9i + = 0, (27) 
Uij - Vij + [\ViVj + Ui k VjVk +i <-> j] = Xij. (28) 



Inserting these into Eq (24), 



dt W ~ Tg t s + T g Q t S g + fXg t p + Vi Q t gi TTij g t Uij, (29) 
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using the notations 

fi^sVi-fP, F^SgVi-Ff, (30) 
aij = (-w + Ts + ViQi + fip + T g Sg)5ij 

+TTij - HikUjk - KjkUik + 9iVj - 0~ij - Efj, (31) 

we obtain 

ViQi = Vi{Tfi + T g Fi + jjji + VjOij - yjiTij) (32) 
-R + f,'\,T + yiVjirij + rfvij + XyTTy + "7; 

+ Sg % + i^V^ - 7 T 9 2 

and deduce 

Qi = Tfi + TgFi + Mi + Vjdij - Ujitij, (33) 
R = fPViT + yiVjTTii 

+aP j v ij +X ij n ij + 'ftf, (34) 
R g = Egv y + f?V,T, - <ylj. (35) 

Given the expressions for R, we may take flux vector as Z = (fP,yi,(jP,Xij), the force 
vectors as Y = (VjT, Vj7Ty, u^-, 71"^), and again formulate the Onsager force-flux relation 
as Z = c ■ F. Analogously, given R g , we have Z g = c g ■ Y g , where Z g = (Ff,T,P) and 
Y g = (ViT g ,Vij). In addition, we require 

X i:j ^ for T g -> 0, (36) 

to ensure permanent elasticity in granular statics. 

This completes the derivation and presentation of the structure of a hydrodynamics of 
permanent elasticity at T g = 0, and transient elasticity at finite T g . To find Granular solid 
hydrodynamics, we still need to specify the energy w, and the functional dependence of the 
transport matrices, c,c g . Instead of a microscopic derivation of these quantities starting 
from some specific interaction, we employ general considerations (such as requiring w to 



have a positive curvature where the system is stable, see § III A) and experimental data to 



narrow down the possibilities. Hereby, w may be determined by static data alone, but c, c„ 
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must be considered using data from granular dynamics. The simplest example is again given 
by c, c g being both diagonal, 

fP = KViT, FP = K g ViT g , y i = P P V j >Ki j , (37) 
= CgViihj + Vgtfj, o% = (vndij + jjv^, (38) 
*a = + PA^u- (39) 



In the next section, § III , an energy expression appropriate for granular media is presented, 
and shown to account for important features of granular statics. For the homogeneous case, 
with VjT, VjT g , VjTtij = 0, we propose to combine this w with the following transport 
structure, diagonal except for the two terms preceded by a, 

erg + Eg = (C + QvuSij + ( V + n g )v% + anj, (40) 

c 

v u ij UuOij . . 

Xij = -av lj 1 J -. (41) 

T Tl 

The first equation is simply a sum of the two dissipative stress contributions. The second 
equation uses the specific form of w, a result of which is 



VA(BA 5 tj - 2Au%) + A—t=6 ij} (42) 



see Eq ( 52 ) below. So the relaxation times are given as 



^2/WA, I=3/WE(b+^|V (43) 

Obviously, a simplification is given by taking either f3 and ft, or r and T\, as independent 
from Uij. Choosing the second possibility, and taking r, T\ as proportional to T g , all other 
coefficients (ie. C^XgiViVgi 01 ) as constant gives us a complete and well specified theory. 
As will be shown in an accompanying paper [5j, this choice leads to a surprisingly good 
agreement with hypoplasticity j6], a modern engineering theory widely employed to model 
solid granular behavior, especially triaxial experiments. 

1. Granular Gas 

Since we are considering a hydrodynamic theory, we should expect the equations as given 
above to easily connect to that of granular gas, such as given in [3] by Haff. Taking the 
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elastic strain to relax infinitely fast, t,t\ — > 0, essentially eliminates Uij as an independent 
variable. As a result, we have w = w(T,T g ,p) in the rest frame, and only Eqs p5|26|27[ ) 
remain as equations of motion, with the dissipative currents given by the second of Eqs (l37|), 



and the first of Eqs (38). Following Haff, we may take w ~ T g , s g ~ hiT g , and the term 



(T g Sg + pp — w) 5ij as the main contribution to the pressure [see Eq (31 )]; also 



(g, %, KgTg, jTg ~ Q ^Tg . (44) 



[Because s is not included as an independent variable, the first of Eqs (26), is ignored in [3], 



as are k, ?7, C Moreover, are included only in R g , not in the stress flux <Tjj, which is 
perhaps not quite consistent. The general gist, however, is certainly the same.] 



III. A GRANULAR ENERGY EXPRESSION 

Linear elasticity is a simple, consistent and complete theory. It starts with an energy w 
that depends on the strain, iiy = |(VjC/j + V jUi), with C/j the displacement vector, 



w 



\Ktf + ^ul (A=-t*a,u.= yk°4), (45) 



see [Mj. K, fi > are two material-dependent constants, referred to as the bulk and shear 
modulus, [uu is the trace of w^-, and tt°- = — |m« 5^ its traceless part.) The stress-strain 
relation is obtained derivative, 

crij = Uij = -^77 = K ^ $ij - u^, (46) 
which contains the pressure P and the scalar shear stress cr s , 



P = \a u = KA, a s = x ja%a% = 2/iu s , (47) 



both employed frequently below. Note that as there is no difference between and tt^ in 



statics, we shall use them interchangeably here, in § [ITT 

Some ramifications of linear elasticity are: (1) Since the stress is given as a function 
of three variables, Ui, the three components of the force balance V^-er^- = pd (with p the 
density and Gi the gravitational constant) suffice to uniquely determine Ui, from which the 
stress <7jj may be calculated for arbitrary geometry. (2) The inverse compliance tensor, M^i, 
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linking the increments of stress and strain, dcr^ and duke, is both isotropic and constant, 

dcr 

d(T *i = "a ~ = Mi M ( 48 ) 
OUkl 

M ijk i = KSijSki - n(5 ik 5ji + 5 jk Sie). (49) 

(3) As the pressure P = KA does not depend on the shear u s , there is no volume dilatancy, 
(dP/du s )\A = 0. (4) Yield is not predicted. [Note that while the points (2), (3), (4) depend 
on the form of the energy w, the statement under (1) is quite general.] 

These equations account well for ordinary solids, but not for granular systems. Sand 
displays volume dilatancy, possesses a compliance tensor with significant stress-induced 
anisotropy, and most importantly, never strays far from yield, displaying significant irre- 
versible, fluid-like, plastic movements in its vicinity. 

The first attempt to modify linear elasticity, so as to better account for granular behavior, 
was due to Boussinesq [15]. He assumed, around 1874, stress-dependent elastic moduli, 



K,fi~ A 1 ' 2 _ pi/3 > in Eq (46) 



nr ( a r 3 — 6u n \ 3 — Qu 2u 
<^VA(A^-— u%y — = ±, (50) 

with v the constant Poisson ratio. This nonlinear stress-strain relation, sometimes referred 
to as the "quasi-elastic model," is employed to understand granular compression [16J and 
sound velocity [17]. Unfortunately, the above failure list of linear elasticity remains partly 
intact: • As P remains a function of A alone, dilatancy vanishes, dP/du s \& = 0. • Yield 



must still be postulated. In addition, Eq (50) contains a basic deficiency: No energy w 



exists such that <Jij = —dw/duij holds, because the associated Maxwell relation is violated, 
duij/dunk ^ dagk/duij. 

We choose the granular elastic energy to be [1] 

w = VA(|SA 2 + ^), (51) 

with A, B > denoting two material constants. The associated stress is 



u 2 



VA(BA 6 tj - 2A u%) + A—^Sij. (52) 



As compared to Eq (50), the only difference is the last term ~ u s /y A. This is, however, 
amazingly useful in accounting for granular behavior. As we shall see, it yields volume 
dilatancy, shear-induced anisotropy, and above all, predicts yield at the Coulomb condition, 

a s /P = ^2A/B. (53) 
14 



In granular materials, there is a regime in which dissipation is insignificant and elastic 
responses dominant: small-amplitude perturbations from given points in the stress space. 
This is shown by Kuwano and Jardine [IS] experimentally, who observed that stress incre- 
ments become reversible if the strain fluctuations are around 10 -4 . It is also corroborated by 
Alonso-Marroquin and Herrmann [13] in molecular-dynamic simulations: Reducing elastic 
strains to 10~ 6 , the irreversible plastic contributions are found around 10~ 14 , implying a line 
as the stress-strain response, rather than the usual ellipse at higher amplitudes. 



This fact is important because it makes a direct verification of Eq (52) possible: Mea- 



sure d<Tjj = (day / duki) duke and diiki independently, and compare the result to = 



daij/duke as calculated from Eq (52). The data in [18] are extensive, comprising of 36 
independent components of My^, all as functions of pressure, shear and the void ratio e. 
Comparing these data to the calculate M^ke is the main result of this section, and represents 



an ambitious test of the energy w, Eq (51): Energy and stress of Eqs ( 51|52 ) depend only 



on two material parameters, A and B, with their ratio fixed by the yield condition, Eq (53). 
Since the Ham river sand used in the experiment has a Coulomb yield angle of around 28°, 
implying £ = B/A = 5/3, only A, a scale factor and a measure of the total hardness, is left 
as an adjustable parameter. Taking A = 5100 Mpa, we find satisfactory agreement with 
their data at all values of pressure and shear, for the void ratio e = 0.66 — except close to 
yield which, due to increased plastic contributions, represents an especially difficult experi- 
mental regime. Because Kuwano and Jardine noticed that e only alters the total hardness, 
by the factor / = (2.17 — e) 2 /(l + e), taking A,B ~ / achieves agreement with respect to 
any other values of e as well. Similar agreement to their data on ballotini (glass beads) was 
achieved by taking A = 4200 Mpa. Therefore, we take 

, , (2.17- ef ^ B 5 

A = AoX 1 3736(1 +V ^ = 3 (M) 

with .4.0 — 5100 and 4200 Mpa being the value of A for e = 0.66, for Ham river sand and 
ballotini, respectively. 

Given this experimental support on the functional dependence of on Uk, we have 



employed Eq (52) to evaluate static stress distributions in silos, sand piles and under point 



loads, not surprisingly with rather satisfactory results, see [2J. Note that Eq (52) does not 
contain any fit parameters: £ = 5/3 is fixed by the yield angle, while Aq, as a scale factor, 
does not enter the stress distribution at all. (Given a solution, one may change the strain 
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by the factor a, and Aq by a -15 , with the stress unchanged and still a solution, provided 
the boundary conditions are the usual ones, either given in terms of stresses or require that 
the displacement vanishes.) 



A. Yield and Energetic Instability 

A thermodynamic energy must be a convex function of state variables to ensure stability 
- this is why compressibility and specific heat are always positive, cf. [20J. Being a quadratic 



function of A and u s , the energy of linear elasticity, Eq (45), is always convex. Conversely, 



the granular energy, Eq (51), is convex if and only if 



{d 2 w/dA 2 ) Us > 0, (d 2 w/du 2 s ) A >0, (55) 
(d 2 w/dAdu s ) 2 < (d 2 w/dA 2 ) Us (d 2 w/du 2 ) A (56) 

hold. (See appendix on some subtleties in this context.) More explicitly, this implies 

u 2 JA 2 < IB I A, (57) 

drawing the boundary for the region of stable strains. Deriving 4P/a s = (A/u s ) x 



(2B/A + u 2 /A 2 ) from Eq (52), and inserting u 2 /A 2 = 2B/A into it, Eq (53), the Drucker- 



Prager version of the Coulomb yield condition (cf. Schofield & Wroth, 1968; Huang, 1983) 
is obtained. The actual Coulomb yield condition, cr s /P = (a/ 18 + 6L 2 sin Lp c ) / (L sin ip c + 3), 
where L = a/3 tan [1 arcsin (a/6 c r ij&jk a ki/ a s)] denotes the Lode parameter, would only re- 



sult if terms ~ u % u °jk u ki are included in Eq (51) 



In a classic paper, Goddard p2] started from Hertz contacts between grains, and con- 
sidered the structure of the energy and stress. He concluded that, if the topology of the 
grain contacts do not change with stress, the energy is a homogeneous function of degree 
5/2 in the strain Uij, of the form w = A 2 - 5 x g(u 2 /A 2 u^u^ k u < l i / A 3 ), where g is an arbitrary 
function. As Eq (7) is clearly a special case of this general energy, we take this as a further, 
microscopically founded support for our starting point. 



There is an instructive analogy between the granular stress-strain relation, Eq (52), and 
the van der Waals equation of state for real gases. The Boyle's law is stable everywhere while 
the van der Waals equation has a non-physical zone, the liquid-gas instability, in which the 
compressibility is negative. Similarly, the Hooke's law is stable everywhere, but the granular 
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a s : Mpa 




FIG. 1: Shear stress versus shear strain for given pressure: for granular elasticity, linear elasticity 
(upper insert), and elastoplastic theory (lower insert). 

stress-strain relation has a forbidden region, that of yield. Note 

dP/dAL > (58) 



is implied by Eqs ( 55|56 ), see appendix, so this forbidden region is also characterized by a 



negative compressibility. The actual innovation of the van der Waals theory is the fact that 
the condition for the onset of the liquid-gas transition, instead of being an extra input, is 
implied by the free energy. Similarly, yield is now a result of elasticity. 

B. Granular Stress-Strain Relation 



The granular stress-strain relation, Eq (52), and the definitions of Eq (47) imply 



P = A 3 l 2 (B+\Au 2 JA 2 ), (59) 
a s = 2AA 1/2 u s . (60) 

Eliminating A, we obtain 

Bo\ - 8A 3 Pu 3 s a s + 8A 5 u 6 s = 0. (61) 
Fig. [1] plots a s versus u s for the fixed pressure of P = 0.1 Mpa. Note how remarkably linear 
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FIG. 2: Thick line: Pressure versus compression at fixed shear. Dashed lines represent unstable 
states. Thin straight line: The same curve for linear elasticity. Insert: The analogous instability 
in the isothermal curve of the van der Waals equation of state. 

the plot is - almost until yield, where the curve turns back abruptly. (Dashed lines are used 
throughout for unstable states.) This behavior is approximated by the elastoplastic model, 
frequently used in soil mechanics: Linear elasticity followed by yield and flat plastic motion, 
see the lower inserts in Fig. [1} Nonlinearity is relevant only when yield is close. 
If instead u s is eliminated from Eqs( 59p0 ), the expression 



o 2 s + 8ABA 3 - 8APA 3 / 2 = (62) 

allows a plot of pressure P versus compression A, at given a s = 0.1 Mpa, see Fig. [2} 
The pressure increases with the compression, implying a positive compressibility, only in 
the region of large A. The compressibility is negative where A is small, and the stability 



condition, Eq (53) or (58), is violated. The van der Waals equation of state, (P — a/v 2 ) (v — 
b) = RT, is quite similar, where \ jv corresponds to A, R is the gas constant and v the 
molar volume, see eg. [2D]- The system can be either in the dense liquid state or the rarefied 
gaseous phase, with the zone in between forbidden, see insert of Fig. [2} 

Alternatively, we may plot A versus u s at fixed P, or P versus a s at fixed A, see Figs. [3] 
and[4j both showing clear evidence of "volume dilatancy," the fact (first noticed by Reynold) 
that granular systems expand with shear, or dA/du s \p ^ 0, or dP/da s \& ^ 0. For linear 
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elasticity, these plots are simply horizontal, and the derivatives vanish. If the Boussinesq 



model, Eq (50 ), were employed, all four plots would be indistinguishable from those of linear 



elasticity. So the last term of Eq (52) is indeed essential. (Plastic motion, not considered 
here, contribute to additional dilatancy, and may dominate.) 



C. Shear-Dependence of the Elastic Moduli 



The Hooke's law, Eq (46), = KA5ij — 2//w°-, may be written as 



2/i' 



with the Poisson ratio v and the Young modulus E given as 

9uK 3K - 2a 

E = — — , v 



(63) 



3K + n' 



6K + 2fi 



(64) 



Requiring the granular stress-strain relation Eq (52) to assume these familiar forms, either 



Eq (46) or (63), leads to strain-dependency of K, fi, 



K = A 1 / 2 (B + lAu 2 JA 2 ) 



(65) 
(66) 




FIG. 3: Compression A versus shear strain u s , at fixed pressure. The dashed line is again unstable. 
In linear elasticity, the same curve is a horizontal straight line. 
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FIG. 4: Pressure P versus shear stress a s , at fixed compression. The dashed line is unstable. In 
linear elasticity, the same curve is a horizontal straight line. 



and via Eq (64) also of E, v. As this is an intuitive way to characterize nonlinear elastic 



behavior, we shall consider their shear and pressure dependency more closely here. Using 
Eqs ( 59p0 ), we write these moduli as 



E 



i4 



3B - 2A£ 



3B + AC 



v 



(67) 



where ^ quantifies shear, 



l±Jl-(B/2A) (a s /PY 



(68) 



and n, K, E, v denote the respective value without shear, at £ = 1, 



ll = A[-\ , K = B 



B 



B 



E 



9AB (P 



3B + AKB 



(69) 



see Fig. [5j (The positive sign in Eq (68) is the stable branch, which meets the unstable 
branch with the negative sign at yield, where the square root vanishes.) 

As mentioned in the introduction, the P 1/,3 -dependence of the twiddled letters is well- 
known. For typical granular behavior, however, the more relevant dependence is that on 



shear, which derives - same as yield and dilatancy - from the last term of Eq (52) 
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yield 



FIG. 5: Variations of K,[i,E,v with cr s /P. The moduli are rescaled by their values at 
denoted respectively with a twiddle. Their variation ~ p 1 / 3 is shown in the insert. 

D. The Compliance Tensor 

1. Theoretical Expressions 



Starting from Eq (52), the tensor Mij k e of Eq (48) is calculated as 



M m = AVA [(ui/4A 2 + 4/3 - ?,B/2A)5 ij 5 kl 
SikSji - Su5 jk + (uij5 k i + 5ijU kl )/A). 



The compliance tensor \j k £, defined via 



duij — Xij k £<ia k £, 
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A,-,* 



(72) 



is obtained by inverting Mijke, 

[Au 2 s + 2(A-B) A 2 } SktSjj _ 6 ik S je + S ie S jk 
Ki3ke ~ 6AA 1 / 2 (Au 2 - 2BA 2 ) AAA 1 / 2 
UjjA&ki + UkjASjj + UjjUkj 
3AV2 ( Au 2 _ 2BA 2 ) ' 

9,4 5 <x 2 + 8 (44 - 9g) ^ 6 _ SjkSjj + 
54// (^a 2 - 8/i 6 i3) w ij 4/i 

4-4 3 /i 3 (44, + alAj) ~ 3^V|4 

9/i (-4V 2 - 8/i 6 £) ' 1 J 

In the first expression A^ is strain-, in the second stress-dependent — where the conversion 

is calculated using A = fi 2 /A 2 , = — |cr^-//x, w s = |cr s //i, with // = ^(^P/iS) 1 / 3 , cf. 

Eqs ( 67p9 ). The second expression - a surprisingly complicated one if the starting expression 



for the energy serves as a benchmark - is what may be compared to experiments directly. 



Before we do this, it is useful to pause and notice that the last term of both Eq (70) 



and (73) deviates structurally from the isotropic form of Eq (49). More generally, for an 



isotropic medium and in the presence of pure compression (a^ = 0, P ^ 0), we may (quite 



independent of the specific form of the elastic energy) take \jkt to be 

^ijki = ^ijhi + ^(SikSje + StfSjk), (74) 

with Ai, A2 arbitrary scalar functions of A, u s , and the Lode parameter L. This is because 
• both Oij and uu are symmetric, hence = Xjiu = \jtk] • the Maxwell relation holds, 
d 2 w /diiijduik = d 2 w/duikduij, hence Ay&« = XkMj- In the presence of shear, 7^ 0, \ijki 
can take on many more terms. To linear order in cr?-, these are 

To second order, we may substitute all above cr°- with er° fc crjL, and also add the terms: 
a ij a k£ anc ^ Vihtfi + a % a u- We shall refer to A°- fc£ as being isotropic, and the cx°-dependent 
ones as displaying "shear-induced anisotropy." If the medium were inherently anisotropic, 
say because the grains are pressed into some quasi-periodic array, leading to a preferred 



direction n, the above expression is more complicated, because 8^ in Eq (74) may now be 
substituted by three different tensors: 5ij — nirij, riiUj, and e^kUk- For triclinic symmetry and 
without the Maxwell relation, all 36 elements of \jki are independent - even in the absence 
of shear. As mentioned, this "fabric anisotropy" is not included in the present consideration, 



because the starting expression for the energy, Eq (51), is isotropic. 
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2. Comparison with Experiments 



Because and symmetric, each characterized by six independent components, 



Eq (71) may be written as a vector equation, du = Ada, with A a 6x6 matrix, and du, da 



given as in Eq (75). In the so-called "principle system" of coordinates, in which is 



diagonal (but not Saij), Kuwano and Jardine take this vector equation to be given as [18 



' dun 
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\ ~da 12 J 



(75) 



with 

f -1/E 1 u 12 /E 2 u 13 /E 3 ^ 

C= V2l /E 1 -I/E2 v 23 jE 3 • (76) 

yVzx/E-i v 32 /E 2 -1/E 3 j 

Gij is referred to as the shear modulus in the i — j plane, Ei the Young modulus along i, 
and Vij the Poisson ratio for "the effect of the 2-strain on j-strain." Identifying these moduli 
with components of the \jki tensor, 

Gij l/4Ajjjj, 
E{ 1/Ajjjj, 

= —\ijjj^jjjj (77) 



(for i y£ j and without summation over i or j), we may employ Eq (73) to obtain 



Gi 3 
Ei = 



Vi 



- G 23 — G\ 2 — /1, 

27/i (A 5 a 2 s - 8fi 6 B) 
9A 5 a 2 s -72^B-Asf 
1 9A 5 a 2 - 72fi 6 B + 2AsiSj 



(78) 
(79) 

(80) 



2 9A 5 a 2 s - 72^B - As) ' 

with /i = A^P/B) 1 / 3 , Si = 3A 2 a® — 4/i 3 , a® = a,i — P, and a,i denoting the three diagonal 
components of in the principle system. Before embarking on a comparison, we shall first 
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establish a few qualitative features from theory: • Without shear, of — > 0, all Ei are equal, 

F^F 27AB ( P Y (811 

Ei -> E sec - —— — — , (81) 



(82) 



2.4 + 98 V^, 

where E sec is called the secant Young modulus. Same holds for the Poisson ratios 

_ 1 98 - AA 



v. 



2 98 + 2 A 



(Note v* differs from u, and E sec from i?, by a constant factor.) • Because of Eq (74) 
and irrespective of the energy specified, we have E\ — Ei — E%, v\2 = = z/23, and 
G12 = G13 = G23 in the absence of shear, of^ = 0. Any discrepancy with experiment 
therefore implies fabric anisotropy. • Finite shear will split Ei and Uij, but not Gij, cf. 



Eq (78) — though this is an energy-related feature. • Because of the Maxwell relation, the 



matrix A of Eq (75) is symmetric, implying especially (no summation) 

VijEi = VjiEj. (83) 
This symmetry was noted by Love (1927) and adopted by Kuwano and Jardine in interpret- 



ing their data [18]. • The moduli E } p } iy are related as E = 2/i {y + 1), see Eqs.(64). A 
similar relation holds for /1, E^ v ik [no summation, see Eqs fl79|80| )], 



Ei - E 3 f = AE 3 (3/i - Ei) (3// - Ej) (84) 

It is important to realize that all formulas of this section hold not only for Cartesian 
coordinates, i — > x,y, z, but also for cylindrical ones, i — > z, p, <p. Taking A = u w + 



u pp + u zz , and similarly for u s , we may again start from the same energy, Eq (51), and 
derive all the results here. [Spatial differentiation is what mars the similarity. Yet once 
the strain components u pp , u ptp . . . are given, no spatial differentiation is needed.] The one 
difference is, for any constant Oij in Cartesian coordinates, there is always a principle system. 
In cylindrical coordinates, this holds only if the stress is also cylindrically symmetric. In 
other words, only if the stress is uniaxially diagonal, a i3 = diag ( <j\, cr 2 , cr 3 ) with a 2 = Q\ in 
Cartesian coordinates, will it be diagonal cylindrically. 

Because Kuwano and Jardine [18] used an axialsymmetric device for their measurements, 
the stress they apply is indeed cylindrically symmetric, with: G pz = G vz , E p = E v , v pz = 



Vtpzi v z P = v z<fi , cf. Eqs.(|78|80|) noting s p = s v . In addition, Eq (|83|) leads to v pip = v vp . 



Following them, we refer to the response coefficients being measured as: Ghh = G pv , G v h 
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0.15 



P: Mpa 



FIG. 6: Variation with pressure P of the shear moduli G v h, Ghh, Young moduli E v , Eh and Poisson 
ratios v v h, fhh (insert), at &h/cr v = 0.45. Symbols are the same data on Ham River sand, at a void 
ratio of 0.66 by Kuwano & Jardine, (2002). 



'fi: v -' r "::j Eh — Ep E l p 1 E v — E Zl Vhh — ^pf ^<pp, ^hv — V pz ^fz, ^vh — L'zp ^z-p 



G pz — G 

where h is the horizontal directions, either p or tp, and v the vertical direction z, see the 
cylinder of Fig. [6j The main plots of Fig. M compare the theoretical curve [calculated by 



taking <J p = o v = Oh and a z = o v in Eqs.(78 80)] and the experimental data [measured with 
Ham River sand] of E h , E v , G v h, Ghh, as functions of P, for ah = 0.45 a v . The insert shows 
the same comparison for u V h, v^h- We especially note that theory and experiment agree on 
the ordering of the induced anisotropy, ie v v h > Vhh, E v > Eh and Ghh ~ G v h, which are 
pairwise equal in linear elasticity and the Boussinesq model. (The slight difference between 
Ghh, G v h is, as mentioned, the result of fabric anisotropy present in the sample.) For a theory 
without any useful fit parameter, the agreement must be considered a convincing verification 
of the elastic approach which, instead of postulating the stress-dependence of 21 (or even 36) 
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P: Mpa 



FIG. 7: Variation of Young and shear moduli, E sec and /i, with pressure P, for the case of vanishing 
shear, a v = ah- The dotted lines are the empirical formula of Kuwano & Jardine (2002), for the 
Ham River sand at the void ratio e = 0.66. The split is proof of fabric anisotropy. 

independent components of A^z directly, looks for one appropriate scalar expression for the 
energy w. Even if it is heavy-handedly simplified, a large number of geometric correlation 
is preserved by the mere fact that \jki is obtained via a double differentiation. This must 
be the main reason why the calculated Xijki stands up so surprisingly well when compared 
to the extensive data of [TS] . 

Kuwano and Jardine [IS] employ the following empirical formulas (in Mpa) for the Ham 
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River sand, 

E v = 204/ (a v /P a f 52 (85) 

E h = 174/ K/P a ) a53 (86) 

G vh = 72f(a v /P a f 32 (a h /P a f 2 (87) 

Ghh = 81f(a v /P a )- om (a h /P a ) - 53 (88) 

where P a = 0.1013 Mpa is the atmospheric pressure and / = (2.17 — e) 2 / (1 + e). (/ = 1.3736 
for the void ratio e = 0.66.) Fig. [7]shows the theoretical and experimental values for Eh, E v , 
G v h and Ghh, as functions of P for the isotropic case ah = o v . The fact that E h , E v and G v h, 
Ghh are pairwise different, indicates (as discussed above) fabric anisotropy. Moreover, the 
theoretical curves are ~ P 1 / 3 , yet experimental ones seem to back a larger power: ~ P 1 / 2 . 
As discussed, • this is a known contradiction between Hertz contact and sound data, with 
possible explanations provided by Goddard[T7j and de Gennes [21], • and a question of 
simplicity versus accuracy in the present approach. 

FigJS] displays the effect of shear on different moduli, with ah ^ cr v . The upper, middle 
and lower figures respectively plot the Young moduli E i: the shear modulus \i (both scaled 
by their isotropic values, P S ec ; P), and the Poisson ratios v^. In agreement with the empirical 



formulas Eqs (85 88), E v increases with a s /P, while Eh decreases, in the region away from 
yield. As yield is approached, both drop quickly to zero. This critical, pre-yield behavior is 
clearly absent for the empirical formulas and is of interests for future experiments. In theory, 
G v h, Ghh are equal, decreasing with a s /P moderately, by less than 20%. In experiments, the 
shear moduli are split, with one increasing, the other decreasing. The discrepancy between 
the theory and experiment is in the range from a s /P = to 0.6 within 20%. This need not 
be a result from fabric anisotropy, as a more complicated energy expression will also do. 



Variation of the Poisson ratios v v h, Vhv, Vhh is given by Eq (80). As depicted, h> v h and Uh v 
increase, while Vhh decreases, with a s /P, all being divergent at yield. No empirical formulae 
for the ratios are given in [T5], and the two circles in the plot simply depict the values from 



the insert of Fig [I] However, v h h = E h / (2Ghh) — 1 was assumed to hold by the authors, and 



interestingly, it may be derived by taking i = h, j = h in Eq (84), yielding Vhh = Eh/(2fi) — 1. 



Assuming that both coefficients A, B of Eq (51) are proportional to / of Eqs (85 88), 
agreement between experiment and theory is extended to all values of the void ratio. Com- 
parison was also made to Kuwano and Jardine's data gained using glass ballotini [18J. Taking 
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FIG. 8: Upper, middle, and lower figures show the Young moduli, shear moduli and Poisson ratios 
as functions of cr s /P. The dotted lines present the empirical formulas of Kuwano & Jardine (2002) 
for the Ram River sand , at the void ratio e = 0.66. 

A = 4200, B = |^4. = 7000, we find similar agreement. 
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E. The Elastic Part of Flow Rules 



The increment relation, Eq (48), may also be written in the matrix form da = Mdu, with 



M a symmetric 6x6 matrix, and da, du still given as in Eq (75 ). The determinant, det M = 
9A 5 (2BA 2 - Au 2 s ) A, calculated from Eq j70|, vanishes at the yield surface, Au 2 = 28 A 2 , 



because an Eigenvalue, call it mi, also does. (This is not a coincidence as M is the Jacobian 
matrix of the energy function, which is positive only in the stable region. It may be of 
interest to note that the determinant of the Bousinesq model, det M = 9 A 5 (3B + 44.) A 3 , 
never vanishes.) The associated Eigenvector fh\ points along the direction at which a finite 
deformation du ^ may take place under constant stress dcr^ = 0. We refer to rh% as 
the elastic flow direction, since miUdw is only the elastic contribution of the strain. Setting 
da,- 



, %3 — in Eq (48) and using Au\ = 2£>A , we obtain 



dui 



A )<iA 




B o£ 



24. a, 



' dA. 



The calculated duij — > du is the Eigenvector mi. Remarkably, one can rewrite this equation 
as duij/dA = dg/ddij, or 



fh\ || dg/da, with g = \JB/2A(J S — P, 



(89) 



implying that the elastic flow direction is perpendicular to the yield surface, as defined by 
the equation g = 0. If the plastic contribution to the strain field may be neglected for some 
reasons, this property is referred to as the associated flow rule see [22j 123] . 



APPENDIX A: ENERGETIC STABILITY 



In the main text, we considered the convexity of the energy with respect to the variables 
u s and A. Relevant is the convexity with respect to Uij. As the transformation between 
these two sets of variables is nonlinear, we bear the burden of proof that both are equivalent. 

Thermodynamic stability requires the elastic energy to be a convex function of its six 
strain variables, or linear combinations of them. This means that all eigenvalues of the 
Jacobian matrix d 2 w/dX a dXp are positive. We take: X\ = u xy ,X 2 = u xz ,X 3 = u yz , 
Xt = (u X x ~ u zz ) /2, X 5 = (u xx - 2u yy + u zz ) /(2\/3), X 6 = -u xx - u yy - u zz , with Q = 
u 2 = 2 Yla=i F° r an energy of the form w = w(A, Q) = w(X s , Q) and denoting 
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/ = Adw/dQ, a = d 2 w/dA 2 , b = Ad 2 w/dQdA, c = 16d 2 w/dQ 2 , the Jacobian matrix is 

( f + cX 2 cX x X 2 cXiXs 0X4X4 cX x X 5 bX x \ 

cX ± X 2 f + cX 2 cX 2 X 3 cX 2 X 4 cX 2 X 5 

cXxXs cX 2 X 3 / + cX| cX 3 X 4 CX3X5 

-2 



bX 2 
bX 3 

cXxX 4 cX 2 X 4 cX 3 X 4 f + cX 2 cX 4 X 5 6X4 



cX l X 5 
\bX x 



cX 2 X 5 
bX 2 



cX 3 X 5 
bX 3 



cX 4 X 5 
6X4 



/ + cX 2 bX 5 
6X5 a 



with its six eigenvalues given as /i_ 4 = / and 



± 



f + a cQ 1 
2 4 2 




eg 
2 



+ 2b 2 Q. 



(Al) 



They are all positive if, and only if, / > 0, 2af + acQ — b 2 Q > 0, / + a + cQ/2 > 0, or 
equivalently, 

dw 



dQ 



dw d 2 w ^d 2 w 
" >0 " 4 dQ + d^ + 8Q dQ' >0 - 



d 2 w dw ^d 2 w d 2 w „ 



<9 2 u> \ 2 
dQdAJ 



> 0, 



(A2) 
(A3) 



cWdQ ^<9Q 2 «9A 2 

Because w 2 = Q, or 2u s (dw/dQ) = dw/du si 2u s (d 2 w/dAdQ) = d 2 w/dAdu S) Au s Q x 
(d 2 w/dQ 2 ) = u s (d 2 w/du 2 ) — dw/du s , these conditions are equivalent to Eqs ( 55|56 ), or 

(A4) 



dw „ <9 2 w „ 9 2 w 

M >0 ' 9^ >0 ' « >0 ' 



<9 2 W <9 2 U> 

<9A 2 <9u 2 



> 



d 2 w 
duedA 



(A5) 



For the energy of Eq (51), the inequalities (A4) imply A > 0, B > 0, while Eq (A5) gives the 
yield condition (57). Using P = dw/dA and a s = dw/du s we can also write Eqs ( A4|A5 ) 



as 



dP 
a~A 



(dP_ 

\d~A 



>0, 
da s 



da t 
du t 

> 



>0, 



dP 

du s 



(A6) 
(A7) 



The Maxwell relation dP/du s \ A = da s /dA\ u and the thermodynamic identities 
dP/dA\ Us = dP/dA\ as + dP/da s \ A ■ da s /dA\ Us , dP/du s \ A = dP/da s \ A ■ da s /du s \ A , 
imply an alternative stability condition, 



(dP/dA) c 
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> 0. 



(A8) 
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